This notebook runs a longitudinal-functional analysis on simulated data.
library(tidyverse)
library(MASS)
library(Rcpp)
library(splines)
library(LFBayes)
library(plotly)
### Control parameters
set.seed(999)
source("/Users/johnshamshoian/Documents/R_projects/LFBayes/Rfuns/Example/Simulation_funcs.R")
errorvar <- .025
SS <- 20
TT <- 20
t <- seq(from = 0, to = 1, length.out = TT)
s <- seq(from = 0, to = 1, length.out = SS)
n <- 60
tt <- list()
tt[[1]] <- 1:(TT*SS)
tt <- rep(tt, n)
p1 <- 12
p2 <- 12
q1 <- 3
q2 <- 3
Bt <- bs(t, df = p1, intercept = TRUE)
Bs <- bs(s, df = p2, intercept = TRUE)
Bt1 <- bs(t, df = p1, intercept = TRUE)
Bs1 <- bs(s, df = p2, intercept = TRUE)
### Generate key model quantities
H <- GenerateH(q1, q2)
mu1 <- GenerateMu1(s,t)
Lambda <- Loading.Matern(t, p1, q1, Bt)
Gamma <- Loading.Brown.Bridge(s, p2, q2)
Cov <- kronecker(Bs%*%Gamma, Bt%*%Lambda)%*%H%*%t(kronecker(Bs%*%Gamma, Bt%*%Lambda)) + errorvar * diag(SS * TT)
### Generate data from model
x <- mvrnorm(n, mu = as.vector(mu1), Sigma = Cov)
sx <- sd(x)
mx <- mean(x)
x <- (x-mx)/sx
Smooth_scaled_cov <- (Cov - errorvar * diag(SS * TT)) / sx^2
mu <- (mu1 - mx)/sx
y <- lapply(1:n, function(i) x[i,])
missing <- list()
for(ii in 1:n){
missing[[ii]] <- numeric(0)
}
X <- cbind(rep(1, n))
### Visualize a few trajectories
nsub <- 6
small_data <- tibble(id = numeric(),
func_time = numeric(),
long_time = numeric(),
value = numeric())
small_data <- small_data %>%
add_row(id = rep(1:nsub, each = SS * TT),
func_time = rep(rep(t, SS), nsub),
long_time = rep(rep(s, each = TT), nsub),
value = c(t(x[1:nsub,])))
id.labs <- paste("Subject", 1:nsub)
names(id.labs) <- "1":nsub
small_data %>%
filter(long_time %in% c(s[5], s[10], s[15])) %>%
ggplot(aes(func_time, value)) +
geom_line(aes(linetype = factor(long_time))) +
facet_wrap(. ~ id, labeller = labeller(id = id.labs)) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
strip.background = element_rect(
color="black", fill="#FFFFFF")) +
labs(x = "Functional time", y = "Response", linetype = "Longitudinal time") +
scale_linetype_discrete(labels = round(c(s[5], s[10], s[15]), 2))

### MCMC control parameters
iter <- 100 # Number of iterations
burnin <- 10 # Burnin iterations
thin <- 1 # Thinning for each chain
nchain <- 1 # Number of chains
q1s <- 3 # Number of latent factors for functional dimension
q2s <- 3 # Number of latent factors for longitudinal dimension
alpha <- .05 # Type 1 error
neig <- 3 # Number of eigenfunctions for inference
### Processing
MCMC <- mcmcWeakChains(y, missing, X, Bs1, Bt1, q1s, q2s, iter, thin, burnin, nchain, 1)
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
MCMC_eigen <- eigenLFChains(Bs1, Bt1, MCMC, neig, iter, burnin, nchain, s, t, alpha)
### Prepare data for plotting
postmean <- matrix(MCMC_eigen$postmean, nrow = TT, ncol = SS)
fig <- plot_ly(z = ~postmean, x = s, y = t) %>% add_surface()
eigenfunction_tibble <- tibble(time = numeric(), type = character(),
value = numeric(), bound = character(),
number = numeric())
eigenfunction_tibble <- eigenfunction_tibble %>%
add_row(value = c(MCMC_eigen$eigvecFuncmean,
MCMC_eigen$eigvecFunclower,
MCMC_eigen$eigvecFuncupper),
number = rep(rep(1:neig, each = TT), 3),
bound = rep(c("mean", "lower", "upper"), each = neig * TT),
type = "Functional",
time = rep(t, 3 * neig)) %>%
add_row(value = c(MCMC_eigen$eigvecLongmean,
MCMC_eigen$eigvecLonglower,
MCMC_eigen$eigvecLongupper),
number = rep(rep(1:neig, each = SS), 3),
bound = rep(c("mean", "lower", "upper"), each = neig * SS),
type = "Longitudinal",
time = rep(s, 3 * neig))
### Plot mean surface
axx <- list(
title = "Longitudinal time"
)
axy <- list(
title = "Functional time"
)
axz <- list(
title = "Posterior mean"
)
fig <- plot_ly() %>%
add_surface(z = ~postmean, x = s, y = t) %>%
layout(scene = list(xaxis=axx,yaxis=axy,zaxis=axz, aspectratio = list(x = .9, y = .8, z = 1))) %>%
hide_colorbar()
fig
NA
### Plot eigenfunctions
options(repr.plot.width=6, repr.plot.height=3)
number.labs <- paste("Eigenfunction", 1:neig)
names(number.labs) <- c("1":neig)
eigenfunction_tibble %>%
ggplot(aes(time, value)) +
geom_line(aes(linetype = bound)) +
facet_grid(type ~ number, labeller = labeller(number = number.labs)) +
scale_linetype_manual(values=c("dashed", "solid", "dashed"))+
theme_bw() +
theme(legend.position = "none", strip.background = element_rect(
color="black", fill="#FFFFFF"),
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

NA
### Percent variability explained by first few eigenfunctions in functional direction
100 * rowMeans(MCMC_eigen$eigvalFunc)[p1:(p1-neig+1)]
[1] 56.731567 27.272872 7.160911
### Percent variability explained by first few eigenfunctions in longitudinal direction
100 * rowMeans(MCMC_eigen$eigvalLong)[p2:(p2-neig+1)]
[1] 71.550694 16.300779 4.075942
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICBkZl9wcmludDogcGFnZWQKLS0tCgpUaGlzIG5vdGVib29rIHJ1bnMgYSBsb25naXR1ZGluYWwtZnVuY3Rpb25hbCBhbmFseXNpcyBvbiBzaW11bGF0ZWQgZGF0YS4KCmBgYHtyLCBtZXNzYWdlID0gRkFMU0V9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KE1BU1MpCmxpYnJhcnkoUmNwcCkKbGlicmFyeShzcGxpbmVzKQpsaWJyYXJ5KExGQmF5ZXMpCmxpYnJhcnkocGxvdGx5KQpgYGAKCmBgYHtyfQojIyMgQ29udHJvbCBwYXJhbWV0ZXJzCgpzZXQuc2VlZCg5OTkpCnNvdXJjZSgiL1VzZXJzL2pvaG5zaGFtc2hvaWFuL0RvY3VtZW50cy9SX3Byb2plY3RzL0xGQmF5ZXMvUmZ1bnMvRXhhbXBsZS9TaW11bGF0aW9uX2Z1bmNzLlIiKQplcnJvcnZhciA8LSAuMDI1ClNTIDwtIDIwClRUIDwtIDIwCnQgPC0gc2VxKGZyb20gPSAwLCB0byA9IDEsIGxlbmd0aC5vdXQgPSBUVCkKcyA8LSBzZXEoZnJvbSA9IDAsIHRvID0gMSwgbGVuZ3RoLm91dCA9IFNTKQpuIDwtIDYwCnR0IDwtIGxpc3QoKQp0dFtbMV1dIDwtIDE6KFRUKlNTKQp0dCA8LSByZXAodHQsIG4pCnAxIDwtIDEyCnAyIDwtIDEyCnExIDwtIDMKcTIgPC0gMwpCdCA8LSBicyh0LCBkZiA9IHAxLCBpbnRlcmNlcHQgPSBUUlVFKQpCcyA8LSBicyhzLCBkZiA9IHAyLCBpbnRlcmNlcHQgPSBUUlVFKQpCdDEgPC0gYnModCwgZGYgPSBwMSwgaW50ZXJjZXB0ID0gVFJVRSkKQnMxIDwtIGJzKHMsIGRmID0gcDIsIGludGVyY2VwdCA9IFRSVUUpCmBgYAoKYGBge3J9CiMjIyBHZW5lcmF0ZSBrZXkgbW9kZWwgcXVhbnRpdGllcwoKSCA8LSBHZW5lcmF0ZUgocTEsIHEyKQptdTEgPC0gR2VuZXJhdGVNdTEocyx0KQpMYW1iZGEgPC0gTG9hZGluZy5NYXRlcm4odCwgcDEsIHExLCBCdCkKR2FtbWEgPC0gTG9hZGluZy5Ccm93bi5CcmlkZ2UocywgcDIsIHEyKQpDb3YgPC0ga3JvbmVja2VyKEJzJSolR2FtbWEsIEJ0JSolTGFtYmRhKSUqJUglKiUKICB0KGtyb25lY2tlcihCcyUqJUdhbW1hLCBCdCUqJUxhbWJkYSkpICsgZXJyb3J2YXIgKiBkaWFnKFNTICogVFQpCmBgYAoKYGBge3J9CiMjIyBHZW5lcmF0ZSBkYXRhIGZyb20gbW9kZWwKCnggPC0gbXZybm9ybShuLCBtdSAgPSBhcy52ZWN0b3IobXUxKSwgU2lnbWEgPSBDb3YpCnN4IDwtIHNkKHgpCm14IDwtIG1lYW4oeCkKeCA8LSAoeC1teCkvc3gKU21vb3RoX3NjYWxlZF9jb3YgPC0gKENvdiAtIGVycm9ydmFyICogZGlhZyhTUyAqIFRUKSkgLyBzeF4yCm11IDwtIChtdTEgLSBteCkvc3gKeSA8LSBsYXBwbHkoMTpuLCBmdW5jdGlvbihpKSB4W2ksXSkKbWlzc2luZyA8LSBsaXN0KCkKZm9yKGlpIGluIDE6bil7CiAgbWlzc2luZ1tbaWldXSA8LSBudW1lcmljKDApCn0KWCA8LSBjYmluZChyZXAoMSwgbikpCmBgYAoKYGBge3IsIGZpZy53aWR0aD04LCBmaWcuaGVpZ2h0PTR9CiMjIyBWaXN1YWxpemUgYSBmZXcgdHJhamVjdG9yaWVzCgpuc3ViIDwtIDYKc21hbGxfZGF0YSA8LSB0aWJibGUoaWQgPSBudW1lcmljKCksCiAgICAgICAgICAgICAgICAgICAgIGZ1bmNfdGltZSA9IG51bWVyaWMoKSwKICAgICAgICAgICAgICAgICAgICAgbG9uZ190aW1lID0gbnVtZXJpYygpLAogICAgICAgICAgICAgICAgICAgICB2YWx1ZSA9IG51bWVyaWMoKSkKc21hbGxfZGF0YSA8LSBzbWFsbF9kYXRhICU+JSAKICBhZGRfcm93KGlkID0gcmVwKDE6bnN1YiwgZWFjaCA9IFNTICogVFQpLAogICAgICAgICAgZnVuY190aW1lID0gcmVwKHJlcCh0LCBTUyksIG5zdWIpLAogICAgICAgICAgbG9uZ190aW1lID0gcmVwKHJlcChzLCBlYWNoID0gVFQpLCBuc3ViKSwKICAgICAgICAgIHZhbHVlID0gYyh0KHhbMTpuc3ViLF0pKSkKCmlkLmxhYnMgPC0gcGFzdGUoIlN1YmplY3QiLCAxOm5zdWIpCm5hbWVzKGlkLmxhYnMpIDwtICIxIjpuc3ViCnNtYWxsX2RhdGEgJT4lCiAgZmlsdGVyKGxvbmdfdGltZSAlaW4lIGMoc1s1XSwgc1sxMF0sIHNbMTVdKSkgJT4lCiAgZ2dwbG90KGFlcyhmdW5jX3RpbWUsIHZhbHVlKSkgKwogIGdlb21fbGluZShhZXMobGluZXR5cGUgPSBmYWN0b3IobG9uZ190aW1lKSkpICsKICBmYWNldF93cmFwKC4gfiBpZCwgbGFiZWxsZXIgPSBsYWJlbGxlcihpZCA9IGlkLmxhYnMpKSArIAogIHRoZW1lX2J3KCkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIHZqdXN0ID0gMC41LCBoanVzdD0xKSwgCiAgICAgICAgc3RyaXAuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdCgKICAgICAgICBjb2xvcj0iYmxhY2siLCBmaWxsPSIjRkZGRkZGIikpICsgCiAgbGFicyh4ID0gIkZ1bmN0aW9uYWwgdGltZSIsIHkgPSAiUmVzcG9uc2UiLCBsaW5ldHlwZSA9ICJMb25naXR1ZGluYWwgdGltZSIpICsKICBzY2FsZV9saW5ldHlwZV9kaXNjcmV0ZShsYWJlbHMgPSByb3VuZChjKHNbNV0sIHNbMTBdLCBzWzE1XSksIDIpKQpgYGAKYGBge3J9CiMjIyBNQ01DIGNvbnRyb2wgcGFyYW1ldGVycwoKaXRlciA8LSA1MDAwICMgTnVtYmVyIG9mIGl0ZXJhdGlvbnMKYnVybmluIDwtIDEwMDAgIyBCdXJuaW4gaXRlcmF0aW9ucwp0aGluIDwtIDEgIyBUaGlubmluZyBmb3IgZWFjaCBjaGFpbgpuY2hhaW4gPC0gMSAjIE51bWJlciBvZiBjaGFpbnMKcTFzIDwtIDMgIyBOdW1iZXIgb2YgbGF0ZW50IGZhY3RvcnMgZm9yIGZ1bmN0aW9uYWwgZGltZW5zaW9uCnEycyA8LSAzICMgTnVtYmVyIG9mIGxhdGVudCBmYWN0b3JzIGZvciBsb25naXR1ZGluYWwgZGltZW5zaW9uCmFscGhhIDwtIC4wNSAjIFR5cGUgMSBlcnJvcgpuZWlnIDwtIDMgIyBOdW1iZXIgb2YgZWlnZW5mdW5jdGlvbnMgZm9yIGluZmVyZW5jZQpgYGAKCmBgYHtyfQojIyMgUHJvY2Vzc2luZwoKTUNNQyA8LSBtY21jV2Vha0NoYWlucyh5LCBtaXNzaW5nLCBYLCBCczEsIEJ0MSwKICAgICAgICAgICAgICAgICAgICAgICBxMXMsIHEycywgaXRlciwgdGhpbiwgYnVybmluLCBuY2hhaW4pCk1DTUNfZWlnZW4gPC0gZWlnZW5MRkNoYWlucyhCczEsIEJ0MSwgTUNNQywgbmVpZywKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGl0ZXIsIGJ1cm5pbiwgbmNoYWluLCBzLCB0LCBhbHBoYSkKYGBgCgpgYGB7cn0KIyMjIFByZXBhcmUgZGF0YSBmb3IgcGxvdHRpbmcKcG9zdG1lYW4gPC0gbWF0cml4KE1DTUNfZWlnZW4kcG9zdG1lYW4sIG5yb3cgPSBUVCwgbmNvbCA9IFNTKQpmaWcgPC0gcGxvdF9seSh6ID0gfnBvc3RtZWFuLCB4ID0gcywgeSA9IHQpICU+JSBhZGRfc3VyZmFjZSgpCmVpZ2VuZnVuY3Rpb25fdGliYmxlIDwtIHRpYmJsZSh0aW1lID0gbnVtZXJpYygpLCB0eXBlID0gY2hhcmFjdGVyKCksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB2YWx1ZSA9IG51bWVyaWMoKSwgYm91bmQgPSBjaGFyYWN0ZXIoKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG51bWJlciA9IG51bWVyaWMoKSkKZWlnZW5mdW5jdGlvbl90aWJibGUgPC0gZWlnZW5mdW5jdGlvbl90aWJibGUgJT4lCiAgYWRkX3Jvdyh2YWx1ZSA9IGMoTUNNQ19laWdlbiRlaWd2ZWNGdW5jbWVhbiwKICAgICAgICAgICAgICAgICAgICBNQ01DX2VpZ2VuJGVpZ3ZlY0Z1bmNsb3dlciwKICAgICAgICAgICAgICAgICAgICBNQ01DX2VpZ2VuJGVpZ3ZlY0Z1bmN1cHBlciksCiAgICAgICAgICBudW1iZXIgPSByZXAocmVwKDE6bmVpZywgZWFjaCA9IFRUKSwgMyksCiAgICAgICAgICBib3VuZCA9IHJlcChjKCJtZWFuIiwgImxvd2VyIiwgInVwcGVyIiksIGVhY2ggPSBuZWlnICogVFQpLAogICAgICAgICAgdHlwZSA9ICJGdW5jdGlvbmFsIiwKICAgICAgICAgIHRpbWUgPSByZXAodCwgMyAqIG5laWcpKSAlPiUKICBhZGRfcm93KHZhbHVlID0gYyhNQ01DX2VpZ2VuJGVpZ3ZlY0xvbmdtZWFuLAogICAgICAgICAgICAgICAgICAgIE1DTUNfZWlnZW4kZWlndmVjTG9uZ2xvd2VyLAogICAgICAgICAgICAgICAgICAgIE1DTUNfZWlnZW4kZWlndmVjTG9uZ3VwcGVyKSwKICAgICAgICAgIG51bWJlciA9IHJlcChyZXAoMTpuZWlnLCBlYWNoID0gU1MpLCAzKSwKICAgICAgICAgIGJvdW5kID0gcmVwKGMoIm1lYW4iLCAibG93ZXIiLCAidXBwZXIiKSwgZWFjaCA9IG5laWcgKiBTUyksCiAgICAgICAgICB0eXBlID0gIkxvbmdpdHVkaW5hbCIsCiAgICAgICAgICB0aW1lID0gcmVwKHMsIDMgKiBuZWlnKSkKCgpgYGAKCmBgYHtyfQojIyMgUGxvdCBtZWFuIHN1cmZhY2UKCmF4eCA8LSBsaXN0KAogIHRpdGxlID0gIkxvbmdpdHVkaW5hbCB0aW1lIgopCgpheHkgPC0gbGlzdCgKICB0aXRsZSA9ICJGdW5jdGlvbmFsIHRpbWUiCikKCmF4eiA8LSBsaXN0KAogIHRpdGxlID0gIlBvc3RlcmlvciBtZWFuIgopCgpmaWcgPC0gcGxvdF9seSgpICU+JQogIGFkZF9zdXJmYWNlKHogPSB+cG9zdG1lYW4sIHggPSBzLCB5ID0gdCkgJT4lCiAgbGF5b3V0KHNjZW5lID0gbGlzdCh4YXhpcz1heHgseWF4aXM9YXh5LHpheGlzPWF4eiwgYXNwZWN0cmF0aW8gPSBsaXN0KHggPSAuOSwgeSA9IC44LCB6ID0gMSkpKSAlPiUKICBoaWRlX2NvbG9yYmFyKCkKZmlnCgpgYGAKCmBgYHtyfQojIyMgUGxvdCBlaWdlbmZ1bmN0aW9ucwoKb3B0aW9ucyhyZXByLnBsb3Qud2lkdGg9NiwgcmVwci5wbG90LmhlaWdodD0zKQpudW1iZXIubGFicyA8LSBwYXN0ZSgiRWlnZW5mdW5jdGlvbiIsIDE6bmVpZykKbmFtZXMobnVtYmVyLmxhYnMpIDwtIGMoIjEiOm5laWcpCmVpZ2VuZnVuY3Rpb25fdGliYmxlICU+JQogIGdncGxvdChhZXModGltZSwgdmFsdWUpKSArCiAgZ2VvbV9saW5lKGFlcyhsaW5ldHlwZSA9IGJvdW5kKSkgKwogIGZhY2V0X2dyaWQodHlwZSB+IG51bWJlciwgbGFiZWxsZXIgPSBsYWJlbGxlcihudW1iZXIgPSBudW1iZXIubGFicykpICsKICBzY2FsZV9saW5ldHlwZV9tYW51YWwodmFsdWVzPWMoImRhc2hlZCIsICJzb2xpZCIsICJkYXNoZWQiKSkrCiAgdGhlbWVfYncoKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiLCBzdHJpcC5iYWNrZ3JvdW5kID0gZWxlbWVudF9yZWN0KAogICAgICAgIGNvbG9yPSJibGFjayIsIGZpbGw9IiNGRkZGRkYiKSwgCiAgICAgICAgYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgdmp1c3QgPSAwLjUsIGhqdXN0PTEpKQogIApgYGAKCmBgYHtyfQojIyMgUGVyY2VudCB2YXJpYWJpbGl0eSBleHBsYWluZWQgYnkgZmlyc3QgZmV3IGVpZ2VuZnVuY3Rpb25zIGluIGZ1bmN0aW9uYWwgZGlyZWN0aW9uCjEwMCAqIHJvd01lYW5zKE1DTUNfZWlnZW4kZWlndmFsRnVuYylbcDE6KHAxLW5laWcrMSldCmBgYAoKYGBge3J9CiMjIyBQZXJjZW50IHZhcmlhYmlsaXR5IGV4cGxhaW5lZCBieSBmaXJzdCBmZXcgZWlnZW5mdW5jdGlvbnMgaW4gbG9uZ2l0dWRpbmFsIGRpcmVjdGlvbgoxMDAgKiByb3dNZWFucyhNQ01DX2VpZ2VuJGVpZ3ZhbExvbmcpW3AyOihwMi1uZWlnKzEpXQpgYGA=